{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 差分方程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 作n阶差分折线图\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def diff_plot(y, n):\n",
    "    if y.shape[0] < n:\n",
    "        print('Error')\n",
    "        return\n",
    "    for i in range(n):\n",
    "        y = np.array([y[i]-y[i-1] if i != 0 else y[0] for i in range(y.shape[0])])\n",
    "    plt.figure()\n",
    "    plt.ylim(y.mean()-(y.max()-y.min())*3, (y.max()-y.min())*3+y.max())\n",
    "    plt.plot(np.linspace(0, y.shape[0]-1, y.shape[0]), y)\n",
    "    plt.title('{} Order Difference of Data'.format(n))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3wUdf7H8dcnBQIEQgmhSw0dC0SaoCiIgP1EDz07ir2Xs9xZ7vSsp6dnP7uiApYDFWwINgQMIL0k9FCSACH01M/vj5n424sbSLI1w+f5eOSR3ZnvzPe7s7vvnflOE1XFGGOMN8VEugHGGGNCx0LeGGM8zELeGGM8zELeGGM8zELeGGM8zELeGGM8zELehJyIzBSRK8Jc51IRGeI+FhF5Q0TyRGSuO+waEckWkT0i0iScbQsHEakjIp+KSL6ITIp0e0zkWMibColIaxEZLyLbRWSviMwVkdMi3KZ2IqJuOO9xg/ozETnZt5yq9lDVme7TQcDJQGtV7Ssi8cBTwHBVTVTV7eF9FWExGmgGNFHVc8uPFJEHRKRIRHa7f6tE5DkRaVHZCiLx422qzkLe+CUijYEfgUKgB5AMPA28JyKjK5gmLshtONj8GqpqInAU8DXwiYhcWkHZtsA6Vd3rPm8GJABLq9mu2OpMF2ZtgVWqWnyQMhNUtT7QGDgbaA7Mq0rQmxpAVe3vMPkD1gG3A4uAfGACkFBB2b8DS4CYcsP/DKwHxH2uwHVABrDWHXYysMKt4zngO+AKn3lcDiwH8oAvgbY+4343v3L1t3PLxJUbfjuQXdZe97UOA8YCB4ASYA/wPrDXncce4Fu3fFecH4sdwErgPJ95vwm8CEx1px0G1AaeBDa49b4E1HHLDwGygNuAHGALcJnP/OoA/3SXYz7Oj2nZtP2BWcBOYCEw5CDvZzdgplt2KXCGO/xBnB/nIvc1jvUz7QPAu+WGxbp1Puk+bwR8BuS679VnOFtDAA+7y/SAW8dz7vBngI3ALmAeMDjSn/vD/S/iDbC/ML7ZTvDNBVrirL0tB66uoOxs4EE/w9u7AdnFfa5uODZ2wyvZ/YKPBuKBW4Bi3JAHzgIy3YCKA/4CzPKZ///Mz0/97fAf8h3c4d18Xusw9/GlwI8VzQOo5wbTZW6begPbgB7u+DfdMD4OZ+s3AfgXMMVtZ33gU+ARt/wQ9zX/zV0Go4B9QCN3/PM44dzKDdaBOD8arYDtbvkYnB/L7UBTP8sh3l2O9wC1gJOA3T7vywOUC/Fy0/sd77Z5jvu4CXAOUNd9jZOA//qUnYnPj7c77EJ3ujicH7mtVLAiYX/h+bPumsPPs6q6WVV34ATT0RWUS8ZZAy1vi8/4Mo+o6g5V3Y8TUMtU9UNVLcIJw60+Za9yyy9XpyvhH8DRItK2gvlV1mb3f+MqTFPmNJzunDdUtVhV5wMf4fxQlZmsqj+pailQAFwJ3OK2c7f7Osb4lC8C/qaqRao6FWdtt4uIxOBsydykqptUtURVZ6lqAU5ATlXVqapaqqpfA+k4y7S8/kAi8KiqFqrqtzhr2udX4/X72oy7DFV1u6p+pKr73Nf4MHDCwSZW1Xfd6YpV9Z84P15dAmyTCUBQ+1BNjeAbuPtw1ur92Qb465tt4TO+zEafxy19n6uqiojv+LbAMyLyT59hgrMWu97P/Cqrlft/RzWmbQv0E5GdPsPigHd8nvu2qSnO2u08ESkbJjhr5WW26//2h+/DCeVknC2B1RW041wROd1nWDwww0/ZlsBG90enzHr+fzlUVyvcZSgidXH2w4zA6boBqC8isapa4m9iEbkNuMJtnwIN+N8VAhNmFvKmIt8A54jIg+WC5DycwFvlM8z3UqZbgDZlT8RJwTY+4zcCD6vq+IPUXZ1Lo56N0/+9shrTbgS+U9WTD1LGt03bgP043TmbqljXNpx+7I44/d/l2/GOql5ZiflsBtqISIzP+3ME//u+VIm7lXE6znsPTndLF6Cfqm4VkaOBBTg/aFDufRKRwTj7bIYCS1W1VETyfMqbCLDuGlORp3HWwl4TkeYikiAi5wP3AneoakVB/DnQQ0T+4B4dcyPOURtlXgLuFpEeACKSJCK/O8SvskSkmYhcD9wP3F3uB6myPgM6i8hFIhLv/h0rIt38FXbr+A/wtIikuO1oJSKnHKoid9rXgadEpKWIxIrIABGpDbwLnC4ip7jDE0RkiIi09jOrOTg7ge902zsEJ6A/qOqLd6fvhrNTujnO4aXg9MPvB3a6R1vdX27SbJx9IfiUL8bZURsnIvfhfIZMBFnIG7/UOXZ8EE7XwjKcHYC3Ahep6oSDTLcNOBd41J0mFfjJZ/wnwGPAByKyC+cInpHVaOJOEdkLLMbpsz5XVV+vxnxw+5uH4/Spb8bp0noMpz+5In/G2fE5230d31D5vufb3Xb/gtM18hjOUUEbgTNxdqbm4qzZ34Gf76mqFgJn4Cy7bcALwMWquqKSbQD4o4jswTk6ZwrO+9VHVcv2b/wLZ2f6Npwd8V+Um/4ZYLR7ktmzOEdKTcPZmliPs8VSna43E0RS8QqZMcaYms7W5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsOi6jj55ORkbdeuXaSbYYwxNcq8efO2qWpTf+OiKuTbtWtHenp6pJthjDE1ioisr2icddcYY4yHWcgbY4yHWcgbY4yHWcgbY4yHWcgbY4yHWcgbY4yHWcgbY4yHBSXkReQWEVkqIktE5H33OtjtRWSOiGSIyAQRqRWMuowxxlRewCEvIq1wbgyRpqo9cW6BNgbnGtlPq2oqzp3exwZalzHGmKoJVndNHFDHvRNQXZxbwJ0EfOiOfws4K0h1GWOMqaSAQ969x+WTwAaccM8H5gE7fW5knEUFNxgWkXEiki4i6bm5uYE2xxhjjI9gdNc0wrllWXucO7TXw//t3PzegkpVX1HVNFVNa9rU7/V1jDHGVFMwumuGAWtVNVdVi4CPgYFAQ7f7BqA1zr0zjTHGhFEwQn4D0F9E6oqIAENxbvw8AxjtlrkEmByEuowxxlRBMPrk5+DsYJ2Pcwf6GOAVnLvZ3yoimUAT4LVA6zLGGFM1QbmevKreD9xfbvAaoG8w5m+MMaZ67IxXY4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxMAt5Y4zxsKCEvIg0FJEPRWSFiCwXkQEi0lhEvhaRDPd/o2DUZYwxpvKCtSb/DPCFqnYFjgKWA3cB01U1FZjuPjfGGBNGAYe8iDQAjse9UbeqFqrqTuBM4C232FvAWYHWZYwxpmqCsSbfAcgF3hCRBSLyqojUA5qp6hYA93+Kv4lFZJyIpItIem5ubhCaY4wxpkwwQj4O6A28qKrHAHupQteMqr6iqmmqmta0adMgNMcYY0yZYIR8FpClqnPc5x/ihH62iLQAcP/nBKEuY4wxVRBwyKvqVmCjiHRxBw0FlgFTgEvcYZcAkwOtyxhjTNXEBWk+NwDjRaQWsAa4DOcHZKKIjAU2AOcGqS5jjDGVFJSQV9VfgTQ/o4YGY/7GGGOqx854NcYYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYD7OQN8YYDwtayItIrIgsEJHP3OftRWSOiGSIyAT31oDGGGPCKJhr8jcBy32ePwY8raqpQB4wNoh1GWOMqYSghLyItAZOBV51nwtwEvChW+Qt4Kxg1GWMMabygrUm/y/gTqDUfd4E2Kmqxe7zLKCVvwlFZJyIpItIem5ubpCaY4wxBoIQ8iJyGpCjqvN8B/spqv6mV9VXVDVNVdOaNm0aaHOMMcb4iAvCPI4DzhCRUUAC0ABnzb6hiMS5a/Otgc1BqMsYY0wVBLwmr6p3q2prVW0HjAG+VdU/ATOA0W6xS4DJgdZljDGmakJ5nPyfgVtFJBOnj/61ENZljDHGj2B01/xGVWcCM93Ha4C+wZy/McaYqrEzXo0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMs5I0xxsMCDnkRaSMiM0RkuYgsFZGb3OGNReRrEclw/zcKvLnGGGOqIhhr8sXAbaraDegPXCci3YG7gOmqmgpMd58bY4wJo4BDXlW3qOp89/FuYDnQCjgTeMst9hZwVqB1GWOMqZqg9smLSDvgGGAO0ExVt4DzQwCkVDDNOBFJF5H03NzcYDbHGGMOe0ELeRFJBD4CblbVXZWdTlVfUdU0VU1r2rRpsJpjjDGGIIW8iMTjBPx4Vf3YHZwtIi3c8S2AnGDUZYwxpvKCcXSNAK8By1X1KZ9RU4BL3MeXAJMDrcsYY0zVxAVhHscBFwGLReRXd9g9wKPARBEZC2wAzg1CXcYYY6og4JBX1R8BqWD00EDnb4wxpvrsjFdjjPEwC3ljjPEwC3ljjPEwC3ljjPEwT4T83oJiJv+6ieKS0kg3xRhjooonQv7ThZu56YNfGfbUd0xM30iRhb0xxgAeCfnz0trw8kV9qFc7jjs/XMSJT87kvTkbKCy2sDfGHN5EVSPdht+kpaVpenp6tadXVWaszOGZ6Zks3LiTlkkJXD2kI+eltSEhPjaILf1/W/MP8OXSrSzbvIu7R3WlYd1aIanHmFDJ3nWAf3+bQcM6tejXoTF92jaibq1gnCdZM+TsOsC/v80EILVZIp1SEklNqU9yYi2cE/qjn4jMU9U0v+O8FPJlVJUfMrbx7PQM0tfnkVK/Nled0JEL+h5BnVqBh31W3j6+WLKVaUu2Mm993m/DLzuuHfef3iPg+dck+fuKyMjZTfauAob3aEZ8rCc2Dg8bM1bkcNukhewpKKakVCkpVeJihF6tk+jbvjH92zchrV0j6ifER7qpQaeqTJqXxUOfLeNAcSm1Y2PYXVD82/ikOvGkpiS6wV+f1BTnB6BFUkLUhf9hF/JlVJWf12zn2ekZzF6zg+TEWlw5uAMX9m9LvdpVW1NZv30v05ZsZdriLSzMygege4sGjOzZnJG9mvPaj2uZlJ7F17eeQPvkekF7DdFAVdm+t5CM7D1k5u4hM3s3GTl7yMjZQ+7ugt/K3TKsMzcNS41gS6tvX2ExB4qq370XHys1KggLi0t54ssV/OeHtXRr0YDnLjiGZg0SmLc+jzlrtjN37Q4WZu2kqESJEejR0gn9fu0b07d94xq/xZqVt4+7P17MDxnb6NuuMY+e04v2yfXI2V1ARvYeMnKcz3im+zhvX9Fv0ybWjqNjSqLzA+D+CKSm1KdN47oRez2Hbcj7mrt2B//+NoMfMrbRqG48VwzuwMUD2h70i5mZs4cvlmxh6uKtLNviXD35qNZJjOjZgpE9m9POJ8xzdh9gyBMzOT61KS9d1CckryHUVJXsXQVk5Owm0w3xij7knVLKNmudD/nEX7L4dmUOX918/P8sl2hXUFzCc99m8uLM1RSXVv+7EBsjPH7OkZzTp3UQWxca67fv5Yb3F7AoK5+LB7TlnlHd/HZn7i8sYcGGPOas3cGctdtZsGEnBe5+rq7N69OvfWP6dWhC3/aNSU6sHe6XUS2lpcq7c9bz2LQVKHDXyK5c2K8tMTEHXzPfvqfgtxWbzOzdZObuISN7Dzk+Kzm3ndyZG4ZGZiXHQt7H/A15/Ht6BjNW5tIgIY7LB7XnsoHtSaobj6qyMns30xZvZdqSLazK3gNAn7aNGNmzOSN6Nqd1o4p/rZ+dnsFTX69i0tUDOLZd45C+jmApLinl00WbeW/OBlZs2f27zdXO5TZVU5sl0rzB7zdXs3cdYOg/v6N320a8ddmxUbc568/8DXn8+cNFZOTs4exjWnF0m4bVntcnCzaxOncP39x6As0aJASxlcE1ZeFm7vl4MTECj48+ihE9m1d62oLiEhZuzGfu2u3MWbuD9HV57C8qASA1JZFHzzmSPm2j91bOa3L3cNdHi5m7bgeDU5N55A+9Dvp9roz8fUVk5u7mxZmr+SFjG9/ePoRWDesEqcWVZyHvx+KsfJ79NoOvl2VTv3Ycp/Rszvz1eazZthcR6NuusRvsLWieVLkv7b7CYk58cibNk+rwyTUDD7l2EElFJaV8smATz8/IZP32faSmJDKgYxM3zOvTKSWxyjue3vhpLQ9+uoznL+jNqUe2CGHrA7O/sIR/frWS135aS4sGCfzjD70Y0sXvjcsqbe22vYz41/cM6dKUly/y+12LqH2FxTw4ZRkT0jfSp20jnhlzdMABV1RSypJN+cxZu4O3Zq2jbq1Yvrj5+KjbL1NcUsprP67lqa9XUTsuhr+c1p1z+7QO6orIpp37OenJmYzq1YKn/3h00OZbWRbyB7Fs8y6em5HBN8tz6NuuMSN6Nmd4j2ak1K/e2tik9I3c8eEinj3/GM44qmWQWxu4guISPpq3iRdmZpKVt5+erRpww0mpnNytWcA/SsUlpZzx3E9s31vA9NuGkFjF/R7h8PPq7dz18SLWb9/Hhf2P4M8jugatL/3Fmat57IsVvHRhb0b0jJ4fuRVbd3H9ewtYnbuHa4d05JZhnYkLchB/syybK95O5y+nduOKwR2COu9ArNi6izs/XMSirHyGd2/GQ2f1JCVEW1qPf7GCF2au5tPrB9GrdVJI6qiIhXwlqGpQftlLSpXT/v0ju/YXMf22E0J26GZVHSgqYcIvG3npu9VsyT/A0W0acuPQTpzYJSWoazQLNuTxhxdncflx7fnrad2DNt9A7T5QxCPTVvDenA20a1KXR885kv4dmgS1juKSUs58/idydhfwzS0nkFQ3sjtiVZX35m7gb58uo35CPP/649EMSk0OWV2XvvEL89fnMeOOIRHvoy8sLuWFmZk8PyOTBgnxPHhmD07t1SKk3Yi7DxQx5ImZdEpJ5INx/cPaZXmwkI+u7aoICtYbEhsj/OXUbmzauZ+3Zq0LyjwDsb+whFd/WMPxj8/g/ilLadWwDm9f3pdPrh3ISV2bBf2DeMwRjTi/7xG8OWsdyzZX+la/ITVjRQ7Dn/6eD+Zu4MrB7Zl20/FBD3iAuNgYHjvnSHbsLeQfU5cHff5Vkb+/iOvem8+9nyyhX4cmTLtpcMgCHpzvz19P687+ohKe/HJlyOqpjEVZOznjuR/51zcZjOrVgq9vPYHTjmwZ8tCtnxDPzSd3Zs7aHXyzPHrudmohHwLHdUrmxC5NeW5GJjv2FkakDXsLinn5u9UMfvxbHvp8OR2a1uO9K/sx6eoBHN+5aUg/8Hee0oWGdeL5y38XUxrAESuB2rmvkFsn/Mplb/5CYu04PrpmIPee2j0o50pUpGerJK4c3IEJ6RuZlbktZPUczPwNeYx65ge+WprN3SO78ualx9K0fujXrDulJHLpwHZMSN/IYvcw43A6UFTCI9OWc9bzP5G3r5BXL07jmTHH0Lhe+A73HHNsGzo2rccjU5dHzeVVQh7yIjJCRFaKSKaI3BXq+qLFPaO6sbegmGenZ4S13l0Hinju2wwGPfYtj0xbQbcWDZh41QA+GDeAgR2Tw7IJ2bBuLe4e1Y35G3YyMX1jyOvzZ9riLQx76numLNzMjSd14rMbB3HMEeE58uPmYam0T67HXR8vZn9hSVjqBOfwwJe+W815L/2MCEy8egBXndAxrAcA3DgslSb1avHAp0sJZ1fwL+t2MPKZH3j5uzWcl9aGr245gWHdm4Wt/jLxsTHcPbIba7bt5f25G8Jevz8hDXkRiQWeB0YC3YHzRSR6OmpDKLVZfcb0PYJ3Z69nTe6ekNeXv6+Ip79exaBHv+XJr1ZxdJuGfHztQN4Z24++7cN/OOc5vVs5J5l8sSKsWzM5uw9wzbvzuGb8fJon1WbK9YO4dXgXaseFb99IQnwsj/yhFxt27OPpb1aFpc7c3QVc8sZcHp22guE9mvH5jYPpHaYfNV8NEuK585SuzFufx+RfN4elznnrd3D+K7MpKill/BX9ePScI0mqE7n9IUO7pTCgQxP+9U0Guw4UHXqCEAv1mnxfIFNV16hqIfABcGaI64watwzrTO24GB6dtiKk9cxavY1Bj33LM9Mz6NehCZ9eP4g3LusbkS95GRHhobN7sudAMY9OC33/tKry0bwsTn7qe6avyOHOEV3477XH0b1lg5DX7U//Dk04v28bXv1hDYuydoa0rhVbdzHq2R+Yu3YHD5/dk+cv6B3RkBvdpzVHtk7ikWnL2etz3kUobN9TwHXjF9CyYR0+v2Ewx3UK3X6HyhIR7j21G3n7CnlhxupINyfkId8K8N1ez3KH/UZExolIuoik5+bmhrg54dW0fm2uGdKRr5ZlM2fN9pDUsWBDHle8lU7zpASm3jiY/1ycFvbDtyrSuVl9xg5uz8T0LNLX7QhZPTv3FXL5m79w26SFdEpJZOqNg7l2SKegHyZYVXeN7EZyYm3u/HBRyPpnl27O5/xXZhMj8N/rjuNP/dpG/ES0mBjh/tN7kL2rgBdmZoasnpJS5eYJv7JjXyEv/Kl3xI9m8tWzVRJnH9OK139aS1bevoi2JdTfAn+ftv/pqFPVV1Q1TVXTmjZtGuLmhN/YQR1o3iCBf0xdHvSdkMu37OLSN34hObE246/oF7G11oO5aWgqrRrW4d5PloQk6Dbt3M/ol37mp8zt3HdadyZeNYBOKYlBr6c6kurE8/ezerJi625e+X5N0Oe/OCufC/4zhzrxsUy8agDdWkTP+9+nbSPOPqYV//lhLRu2hybkyi5T8rczetCzVXSs2Pi6fXgXBHgiwkcbhTrks4A2Ps9bA+HpqIsSdWrFcscpXViYlc+ni4L30tdu28tFr82lTnws46/oF7ITPAJVt1Yc95/enZXZu3nzp3VBnffKrbs554VZZO86wNtj+3L5oPbERtlZxqf0aM6oXs15ZnoGq4O4b2bBhjwueHU29RPimHDVANo2ib7rBd01sitxMcJDny8L+ry/X5XLM9MzOKd3a/54bJtDTxABLRvW4crBHZj862YWbgxtl93BhDrkfwFSRaS9iNQCxgBTQlxn1Dn7mFb0aNmAx79YyYGiwI+22LRzPxe+OodSVd69om9Er35XGSd3b8bQrik8/c0qNu/cH5R5zlmzndEvzUJRJl09ICTHvQfLA2f0oE58LHd9tCgoW3Pz1u/gotfm0qhuLSZcNSBq3/9mDRK4/qROfLUsmx8ygtcVu3nnfm76YAFdmtXnobN6Rrx76mCuHtKR5MRaPPz58rAebeQrpCGvqsXA9cCXwHJgoqouDWWd0SgmRrh3lHOC1BsBrs3m7i7golfnsGt/EW9f3pdOKfWD08gQEhEeOKMHpar87dPA1+qmLd7CRa/PJaV+bT6+9ji6No+ebgp/UuoncO+p3fhlXR7jAzysbu7aHVz82lya1q/NhKv6R+RiWFUxdlB72japy4OfLgtKd11hcSnXvzefwuJSnv9T75Ce8xAMibXjuOXkzsxdt4OvlmVHpA0h3zOlqlNVtbOqdlTVh0NdX7Qa2CmZoV1TeGFGJtv3FBx6Aj/y9xVx8etz2Zy/n9cvOzYq+yEr0qZxXW44KZUvlm5lxorqnw34zs/ruPa9+fRs2YAPrx4Y9SFX5tw+rRnUKZnHpq1gS371tmZmrd7GJa/PpXlSAh+M60+LpOh/7bXjYvnLqd3JzNnDOz+vD3h+j05bwfwNO3ls9JF0bBod+14O5Y9pbeiUksij01ZE5AQpO+M1jO4e1ZV9RSU8U40TpPYWFHPZm3PJzNnNyxel1ZhLGfu6cnAHOjatx31TllS520pVefLLlfx18lKGdk1h/BX9aRTGMxkDJSL84+xelJQqf/lkSZU33X/M2Mblb/5Cm8Z1+GDcgKi+nHF5w7qlMDg1mae/WVXtFRyAqYu38PpPa7l0YDtOOzL6Lv5XkbjYGO4Z1ZW12/YyfnbgP3RVZSEfRp1S6nN+3zaMn7OhSjvhDhSVcNU78/h1406eHXMMJ3SumUch1YqL4e9n9WTjjv08P6Pyh9YVl5Ty548W8dyMTMYc24aXLuwT9Zvp/hzRpC63De/M9BU5fLpoS6Wnm7kyh8vf+oV2Terx/pX9w3KJgmASEe4/vTv7C0t48qvqnRy2dtte7vxwEUe3acg9o7oFuYWhd2KXFI7r1IRnpmeQvz+8J0hZyIfZzcM6Uyc+ttInSBWXlHLj+wv4MXMbj48+ipG9oucSttUxsGMyZx/Tipe+W12pH7r9hc4P3MT0LG4cmsojf+gV8ePfA3HZce05qnUSD05ZSl4lzgSevjybcW/PIzUlkfev7E+TGnIHpvI6pdTnkoHt+OCXDSzZVLXr2uwvLOGad+cRHys8/6fe1Iqree+/iHDPqG7s3F/EC1VYwQmGmre0arjkROcEqa+XZTP7ECdIlZYqd364iK+WZfPA6d0ZXQNuLVcZZbebu2/ywbstduwt5IJXZzNjZQ4PndWTW0/uHNVHUlRGbIzw6DlHkr+/iL8f4tDCL5du5ep359G1RX3eq2HdU/7cODSVxnVr8cCUql3X5r7JS1iZvZun/3h0jdkH40+Plkmc07s1b/y0jo07wneClIV8BIwd1J6WSQk8/HnFJ0ipKvdPWcrHCzZx28mdufS49mFuZeg0rV+bO0/pwk+Z25my0P+5Axt37GP0S7NYunkXL/ypDxf2bxvmVoZOtxYNuGZIRz6ev4nvVvk/tHDq4i1cN34+PVom8c7YflF1Nmd1JdWJ545TupC+Pq/C9728ib9sZNK8LG44sVPAd++KBrcP70JMDDwexhOkLOQjICE+ljtGdGHxpnwmL9zkt8wTX67kndnrGXd8B64/qVOYWxh6F/Rry5Gtk3jo8+W/u4jTss27OOfFWWzbXcD4K/pV6T6kNcX1J3WiY9N63PPx4t9d32XKws3c8P4Cjm7TkHfG9o3odWiC7dy0No1nNwgAAA2dSURBVPRs1YBHpq5gX+HBr2uzdHM+f528hOM6NeGmYZ3D1MLQap6UwLjBHfh04WYWbMgLS50W8hFy5lGt6NUqiSf8nCD14szVvDBzNef3PYK7R3at8V0U/sTGCA+f1Yttewp4ymdn3KzV2/jjyz8TI8KkqwfWyKOIKqN2XCyPnXMkm/P38+RX/79W98mCLG7+YAF92jbircv7Bu3WhNEiNkZ44PQebN11gBdnVnzxrl0Hirh2/Hwa1o3nmTHHRN2ZzIEYd0JHkhNrh+0EKQv5CImJcXbEbM4/wGs/rv1t+Duz1/PYFys446iWUX82X6B6tU7iov5tefvndSzZlM9nizZz6eu/0DwpgY+vHUiX5tF/olcg0to15qL+bXlz1jrmb8hjYvpGbp24kP4dmvDmZcdSLwrvkRsMae0ac9bRLXn5+zV++6ZVlTsmLSQrbz/PX9A74rcSDLbE2nHcenJn0tfn8eXSrSGvz0I+ggZ0bMKwbs14ceZqtu0p4JMFWdw3eQnDuqXwz/OO8tTaS0VuG96FxvVqc+Xb6dzw/gKOapPEpKsH0LIG72CrijtHdKVFgwSuemced364iEGdknn90mOpW8ubAV/mrpHdKryuzWs/ruVL965WaR7dkjsvrTWdmzknSBUWh/YEKQv5CLt7VFf2FzmHiN0+aRH92zfhuQt6E1+DDxOsiqQ68fz1tG5syT/Ayd2a8c7YfjSsW7OPIqmKxNpxPHx2L3J3F3Bil6b85+K0qLn5eyg1T0rguhM78eXSbH7yuU1i+rodPDptBaf0aMbYQd452KC8uNgY7h7VjXXb9/FuiE+QkkhdNMeftLQ0TU9Pj3Qzwu6+yUt4++f1HN2mIe9e0Y9Ej26mH8zKrbvplJJ4WGy9+JOZs5u2TeodNj/u4JzkN/zp70mIj2HqjYPZub+I0579kdrxMUy5fpCndjj7o6pc/PpcFm/K57vbTwzoCCoRmaeqaf7GHT6fqCh22/Au3HFKF9687NjDMuABujSvf9gGPDgnCx1OAQ/OUWb3ntqNVdl7ePvn9dz8gc8NQDwe8OCcIHX3yG7k7y/iuRmhuxf04fWpilJJdeK57sROh1U3hTEAw7s3Y3BqMn//fBk/Zm7j72f2oEfLmnPhvUB1b9mAc/u05q1Z60N2cxULeWNMxIgI953WnfjYGEb3ac15adF5A5BQum14F2JjhJe/D839YK1P3hgTcdv2FNCkXi1PHzJ8MHPWbOfI1g2rfeG9g/XJH54dwMaYqOK1Y+Grql8I72xm3TXGGONhAYW8iDwhIitEZJGIfCIiDX3G3S0imSKyUkROCbypxhhjqirQNfmvgZ6qeiSwCrgbQES649y0uwcwAnhBRLx/hocxxkSZgEJeVb9yb9YNMBsou+D5mcAHqlqgqmuBTKBvIHUZY4ypumD2yV8OTHMftwI2+ozLcof9joiME5F0EUnPzfV/bW1jjDHVc8ija0TkG8DfBb3vVdXJbpl7gWJgfNlkfsr7PVZTVV8BXgHnEMpKtNkYY0wlHTLkVXXYwcaLyCXAacBQ/f+D7rMA37MaWgOVuxWMMcaYoAn06JoRwJ+BM1TV95zcKcAYEaktIu2BVGBuIHUZY4ypukBPhnoOqA187Z6pNltVr1bVpSIyEViG041znaqWHGQ+xhhjQiCgkFfVCm8+qqoPAw8HMn9jjDGBsTNejTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGwyzkjTHGw4IS8iJyu4ioiCS7z0VEnhWRTBFZJCK9g1GPMcaYqgk45EWkDXAysMFn8Eicm3enAuOAFwOtxxhjTNUFY03+aeBOQH2GnQm8rY7ZQEMRaRGEuowxxlRBQCEvImcAm1R1YblRrYCNPs+z3GH+5jFORNJFJD03NzeQ5hhjjCkn7lAFROQboLmfUfcC9wDD/U3mZ5j6GYaqvgK8ApCWlua3jDHGmOo5ZMir6jB/w0WkF9AeWCgiAK2B+SLSF2fNvY1P8dbA5oBba4wxpkqq3V2jqotVNUVV26lqO5xg762qW4EpwMXuUTb9gXxV3RKcJhtjjKmsQ67JV9NUYBSQCewDLgtRPcYYYw4iaCHvrs2XPVbgumDN2xhjTPXYGa/GGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhAYe8iNwgIitFZKmIPO4z/G4RyXTHnRJoPcYYY6ouoNv/iciJwJnAkapaICIp7vDuwBigB9AS+EZEOqtqSaANNsYYU3mBrslfAzyqqgUAqprjDj8T+EBVC1R1Lc4NvfsGWJcxxpgqCjTkOwODRWSOiHwnIse6w1sBG33KZbnDfkdExolIuoik5+bmBtgcY4wxvg7ZXSMi3wDN/Yy6152+EdAfOBaYKCIdAPFTXv3NX1VfAV4BSEtL81vGGGNM9Rwy5FV1WEXjROQa4GNVVWCuiJQCyThr7m18irYGNgfYVmOMMVUUaHfNf4GTAESkM1AL2AZMAcaISG0RaQ+kAnMDrMsYY0wVBXR0DfA68LqILAEKgUvctfqlIjIRWAYUA9fZkTXGGBN+AYW8qhYCF1Yw7mHg4UDmb4wxJjB2xqsxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYhbwxxniYOHfriw4ikgusr+bkyTj3l41W0d4+iP42WvsCY+0LTDS3r62qNvU3IqpCPhAikq6qaZFuR0WivX0Q/W209gXG2heYaG9fRay7xhhjPMxC3hhjPMxLIf9KpBtwCNHePoj+Nlr7AmPtC0y0t88vz/TJG2OM+T0vrckbY4wpx0LeGGM8rMaFvIiMEJGVIpIpInf5GV9bRCa44+eISLswtq2NiMwQkeUislREbvJTZoiI5IvIr+7ffeFqn1v/OhFZ7Nad7me8iMiz7vJbJCK9w9i2Lj7L5VcR2SUiN5crE/blJyKvi0iOiCzxGdZYRL4WkQz3f6MKpr3ELZMhIpeEsX1PiMgK9z38REQaVjDtQT8PIWzfAyKyyed9HFXBtAf9voewfRN82rZORH6tYNqQL7+AqWqN+QNigdVAB6AWsBDoXq7MtcBL7uMxwIQwtq8F0Nt9XB9Y5ad9Q4DPIrgM1wHJBxk/CpgGCNAfmBPB93orzkkeEV1+wPFAb2CJz7DHgbvcx3cBj/mZrjGwxv3fyH3cKEztGw7EuY8f89e+ynweQti+B4DbK/EZOOj3PVTtKzf+n8B9kVp+gf7VtDX5vkCmqq5R1ULgA+DMcmXOBN5yH38IDBURCUfjVHWLqs53H+8GlgOtwlF3EJ0JvK2O2UBDEWkRgXYMBVaranXPgA4aVf0e2FFusO/n7C3gLD+TngJ8rao7VDUP+BoYEY72qepXqlrsPp0NtA52vZVVwfKrjMp83wN2sPa52XEe8H6w6w2XmhbyrYCNPs+z+H2I/lbG/ZDnA03C0jofbjfRMcAcP6MHiMhCEZkmIj3C2jBQ4CsRmSci4/yMr8wyDocxVPzFiuTyK9NMVbeA8+MOpPgpEy3L8nKcrTN/DvV5CKXr3e6k1yvo7oqG5TcYyFbVjArGR3L5VUpNC3l/a+TljwGtTJmQEpFE4CPgZlXdVW70fJwuiKOAfwP/DWfbgONUtTcwErhORI4vNz4all8t4Axgkp/RkV5+VRENy/JeoBgYX0GRQ30eQuVFoCNwNLAFp0ukvIgvP+B8Dr4WH6nlV2k1LeSzgDY+z1sDmysqIyJxQBLV21SsFhGJxwn48ar6cfnxqrpLVfe4j6cC8SKSHK72qepm938O8AnOJrGvyizjUBsJzFfV7PIjIr38fGSXdWO5/3P8lInosnR39J4G/EndDuTyKvF5CAlVzVbVElUtBf5TQb2RXn5xwB+ACRWVidTyq4qaFvK/AKki0t5d2xsDTClXZgpQdhTDaODbij7gweb2370GLFfVpyoo07xsH4GI9MV5D7aHqX31RKR+2WOcnXNLyhWbAlzsHmXTH8gv65YIowrXniK5/Mrx/ZxdAkz2U+ZLYLiINHK7I4a7w0JOREYAfwbOUNV9FZSpzOchVO3z3c9zdgX1Vub7HkrDgBWqmuVvZCSXX5VEes9vVf9wjv5YhbPX/V532N9wPswACTib+ZnAXKBDGNs2CGdzchHwq/s3CrgauNotcz2wFOdIgdnAwDC2r4Nb70K3DWXLz7d9AjzvLt/FQFqY39+6OKGd5DMsossP5wdnC1CEs3Y5Fmc/z3Qgw/3f2C2bBrzqM+3l7mcxE7gsjO3LxOnPLvsclh1x1hKYerDPQ5ja9477+VqEE9wtyrfPff6773s42ucOf7Psc+dTNuzLL9A/u6yBMcZ4WE3rrjHGGFMFFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONhFvLGGONh/weK+0CEfXeV0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.87372336 0.19410722 0.69574517]\n",
      "[11.         16.         25.         12.         12.         18.\n",
      " 26.         14.         13.         20.         27.         15.\n",
      " 15.         24.         30.         15.         16.         25.\n",
      " 32.         17.         17.5869272  34.47810922 19.16756978 23.65432962\n",
      " 13.72045164 18.39588258 25.20206786 16.07232878 13.31560489 19.52848113]\n"
     ]
    }
   ],
   "source": [
    "# 例15.6\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "y0 = np.array([11, 16, 25, 12, 12, 18, 26, 14, 13, 20, 27, 15, 15, 24, 30, 15, 16, 25, 32, 17])\n",
    "diff_plot(y0, 1)\n",
    "y = y0[8: 20]\n",
    "x = np.array([y0[4:16], y0[0:12], np.ones(12)]).T\n",
    "z = np.matmul(np.linalg.inv(np.matmul(x.T, x)), np.matmul(x.T, y))\n",
    "print(z)\n",
    "\n",
    "pred_y = y0.copy()\n",
    "for i in range(10):\n",
    "    pred_y = np.r_[pred_y, [z[0] * pred_y[i-4] + z[1] * pred_y[i-8] + z[2]]]\n",
    "    \n",
    "print(pred_y)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
